home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / matrix.lha / Matrix / src / Matrix.ccP < prev    next >
Encoding:
Text File  |  1993-07-07  |  20.0 KB  |  902 lines
  1. pushdef(`index', `ifelse($#, 0, ``index'', ``index($*)'')')dnl
  2. pushdef(`shift', `ifelse($#, 0, ``shift'', ``shift($*)'')')dnl
  3. pushdef(`format', `ifelse($#, 0, ``format'', ``format($*)'')')dnl
  4. pushdef(`index', `ifelse($#, 0, ``index'', ``index($*)'')')dnl
  5. pushdef(`shift', `ifelse($#, 0, ``shift'', ``shift($*)'')')dnl
  6. pushdef(`format', `ifelse($#, 0, ``format'', ``format($*)'')')dnl
  7. #include <<T>.Matrix.h>
  8.  
  9. <T>Matrix <T>Array::operator - () const {
  10.   <T>* newX = new <T> [M*N];
  11.   <T>* t = newX;
  12.   <T>* u =    X;
  13.   <T>* v =    X;
  14.   do
  15.     while (u < v + N)
  16.       *t++ = -(*u++);
  17.   while ((u = v += L) < &X[M*L]);
  18.   return <T>Matrix(M, N, newX);
  19.   }
  20.  
  21. dnl 1$ = {*,/,+,-}
  22. define(arithmetic,
  23. `<T>Matrix operator $1 (const <T&> a, const <T>Array& b) {
  24.   <T>* newx = new <T> [b.m()*b.n()];
  25.   <T>* t = newx;
  26.   <T>* u =  b.x();
  27.   <T>* v =  b.x() + b.n();
  28.   do {
  29.     while (u < v)
  30.       *t++ = a $1 *u++;
  31.     u += b.l() - b.n();
  32.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  33.   return <T>Matrix(b.m(), b.n(), newx);
  34.   }
  35.  
  36. <T>Matrix operator $1 (const <T>Array& a, const <T&> b) {
  37.   <T>* newx = new <T> [a.m()*a.n()];
  38.   <T>* t = newx;
  39.   <T>* u =  a.x();
  40.   <T>* v =  a.x() + a.n();
  41.   do {
  42.     while (u < v)
  43.       *t++ = *u++ $1 b;
  44.     u += a.l() - a.n();
  45.     } while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
  46.   return <T>Matrix(a.m(), a.n(), newx);
  47.   }
  48.  
  49. <T>Matrix operator $1 (const <T>Array& a, const <T>Array& b) {
  50.   if (a.n() == b.n()) {
  51.     if (a.m() == b.m()) {        // a.n() == b.n()
  52.       <T>* newx = new <T> [b.m()*b.n()];
  53.       <T>* s = newx;
  54.       <T>* t = newx;
  55.       <T>* u =  a.x();
  56.       <T>* v =  b.x();
  57.       while ((t += b.n()) <= &newx[b.m()*b.n()]) {
  58.     while (s < t)
  59.       *s++ = *u++ $1 *v++;
  60.     u += a.l() - a.n();
  61.     v += b.l() - b.n();
  62.     };
  63.       return <T>Matrix(b.m(), b.n(), newx);
  64.       };
  65.     if (a.m() == 1) {            // a.n() == b.n() && a.m() != b.m()
  66.       <T>* newx = new <T> [b.m()*b.n()];
  67.       <T>* s = newx;
  68.       <T>* t = newx;
  69.       <T>* v =  b.x();
  70.       while ((t += b.n()) <= &newx[b.m()*b.n()]) {
  71.     <T>* u = a.x();
  72.     while (s < t)
  73.       *s++ = *u++ $1 *v++;
  74.     v += b.l() - b.n();
  75.     };
  76.       return <T>Matrix(b.m(), b.n(), newx);
  77.       };
  78.     if (b.m() == 1) {    // a.n() == b.n() && a.m() != b.m() && a.m() != 1
  79.       <T>* newx = new <T> [a.m()*a.n()];
  80.       <T>* s = newx;
  81.       <T>* t = newx;
  82.       <T>* u =  a.x();
  83.       while ((t += a.n()) <= &newx[a.m()*a.n()]) {
  84.     <T>* v = b.x();
  85.     while (s < t)
  86.       *s++ = *u++ $1 *v++;
  87.     u += a.l() - a.n();
  88.     };
  89.       return <T>Matrix(a.m(), a.n(), newx);
  90.       };
  91.     };
  92.   if (a.n() == 1) {            // b.n() != a.n()
  93.     if (a.m() == b.m()) {        // b.n() != a.n() == 1
  94.       <T>* newx = new <T> [b.m()*b.n()];
  95.       <T>* s = newx;
  96.       <T>* t = newx;
  97.       <T>* u =  a.x();
  98.       <T>* v =  b.x();                        
  99.       while ((t += b.n()) <= &newx[b.m()*b.n()]) {
  100.     while (s < t)
  101.       *s++ = *u $1 *v++;
  102.     u++;
  103.     v += b.l() - b.n();
  104.     };
  105.       return <T>Matrix(b.m(), b.n(), newx);
  106.       };
  107.     if (a.m() == 1) {            // b.n() != a.n() == 1 && a.m() != b.m()
  108.       <T>* newx = new <T> [b.m()*b.n()];
  109.       <T>* s = newx;
  110.       <T>* t = newx;
  111.       <T>* u =  a.x();
  112.       <T>* v =  b.x();
  113.       while ((t += b.n()) <= &newx[b.m()*b.n()]) {
  114.     while (s < t)
  115.       *s++ = *u $1 *v++;
  116.     v += b.l() - b.n();
  117.     };
  118.       return <T>Matrix(b.m(), b.n(), newx);
  119.       };
  120.     };
  121.   if (b.n() == 1) {            // a.n() != b.n() && a.n() != 1
  122.     if (a.m() == b.m()) {        // a.n() != b.n() == 1
  123.       <T>* newx = new <T> [a.m()*a.n()];
  124.       <T>* s = newx;
  125.       <T>* t = newx;
  126.       <T>* u =  a.x();
  127.       <T>* v =  b.x();
  128.       while ((t += a.n()) <= &newx[a.m()*a.n()]) {
  129.     while (s < t)
  130.       *s++ = *u++ $1 *v;
  131.     u += a.l() - a.n();
  132.     v++;
  133.     };
  134.       return <T>Matrix(a.m(), a.n(), newx);
  135.       };
  136.     if (b.m() == 1) {            // a.n() != b.n() == 1 && a.m() != b.m()
  137.       <T>* newx = new <T> [a.m()*a.n()];
  138.       <T>* s = newx;
  139.       <T>* t = newx;
  140.       <T>* u =  a.x();
  141.       <T>* v =  b.x();
  142.       while ((t += a.n()) <= &newx[a.m()*a.n()]) {
  143.     while (s < t)
  144.       *s++ = *u++ $1 *v;
  145.     u += a.l() - a.n();
  146.     };
  147.       return <T>Matrix(a.m(), a.n(), newx);
  148.       };
  149.     };
  150.   a.error("nonconformant <T>Array $1 operands.");
  151.   return <T>Matrix();
  152.   }
  153. ')dnl
  154. arithmetic(*)
  155. arithmetic(/)
  156. <T>Matrix operator % (const <T>Array& a, const <T>Array& b) {    // inner product
  157.   if (a.n() != b.n())
  158.     a.error("nonconformant <T>Array % operands.");
  159.   <T>* newx = new <T> [a.m()*b.m()];
  160.   <T>* s = newx;
  161.   <T>* t =  a.x();
  162.   <T>* u;
  163.   <T>* v;
  164.   do {
  165.     v = b.x();
  166.     while (v < &b.x()[b.m()*b.l()]) {
  167.       u = t; *s = *u++ * *v++;
  168.       while (u < t + a.n())
  169.     *s += *u++ * *v++;
  170.       ++s;
  171.       v += b.l() - b.n();
  172.       };
  173.     } while((t += a.l()) < &a.x()[a.m()*a.l()]);
  174.   return <T>Matrix(a.m(), b.m(), newx);
  175.   }
  176.  
  177. arithmetic(+)
  178. arithmetic(-)
  179. ifelse(basetype, complex, `
  180. extern int Array_number; extern char* Array_format; extern char* Array_space;
  181. ', `
  182. int    Array_number     = 4;
  183. char*    Array_format     = "%g";
  184. char*    Array_space     = " ";
  185. char    empty[1]     = "";
  186.  
  187. char* format(char* s, int n, char* w) {
  188.   Array_number    = n;
  189.   Array_format    = s;
  190.   Array_space    = w;
  191.   return empty;
  192.   }
  193. ')dnl
  194. ostream& operator << (ostream& s, const <T>Array& b) {
  195.   <T>* t = b.x();
  196.   <T>* u = b.x();
  197.   <T>* v = b.x() + b.n();
  198.   do {
  199.     do {
  200.       t += (0 < Array_number)? Array_number: b.n();
  201.       t  = (t < v)? t: v;
  202.       while (u < t)
  203.     if (!(ifelse(basetype, complex,
  204.          `s << "("  && s << form(Array_format, real(*u)) &&
  205.           s << ", " && s << form(Array_format, imag(*u)) &&
  206.           s << ")"',
  207.          `s << form(Array_format, *u)')
  208.           && ((++u < t) ? s << Array_space : s << "\n")))
  209.       return s;
  210.       } while (t < v);
  211.     t = u += b.l() - b.n();
  212.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  213.   return s;
  214.   }
  215.  
  216. istream& operator >> (istream& s, const <T>Array& b) {
  217.   <T>* u = b.x();
  218.   <T>* v = b.x() + b.n();
  219.   do {
  220.     while (u < v && s >> *u)
  221.       u++;
  222.     if (u < v)
  223.       return s;
  224.     u += b.l() - b.n();
  225.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  226.   return s;
  227.   }
  228.  
  229. <T>Matrix operator << (const <T>Array& a, int n) {
  230.   <T>* newx = new <T> [a.m()*a.n()];
  231.   <T>* t =  a.x() + n;
  232.   <T>* u = newx;
  233.   <T>* v = newx;
  234.   while ((v += a.n()) <= &newx[a.m()*a.n()]) {
  235.     while (u < v - n)
  236.       *u++ = *t++;
  237.     while (u < v)
  238.       *u++ = 0.0;
  239.     t += a.l() - a.n() + n;
  240.     };
  241.   return <T>Matrix(a.m(), a.n(), newx);
  242.   }
  243.  
  244. <T>Matrix operator >> (const <T>Array& a, int n) {
  245.   <T>* newx = new <T> [a.m()*a.n()];
  246.   <T>* t =  &a.x()[a.m()*a.l()];
  247.   <T>* u = &newx[a.m()*a.n()];
  248.   <T>* v = &newx[a.m()*a.n()];
  249.   while (newx <= (u -= a.n())) {
  250.     t -= a.l() - a.n() + n;
  251.     while (v > u + n)
  252.       *--v = *--t;
  253.     while (v > u)
  254.       *--v = *t;
  255.     };
  256.   return <T>Matrix(a.m(), a.n(), newx);
  257.   }
  258.  
  259. dnl 1$ = {==,< }
  260. define(comparison,
  261. `int operator $1 (const <T>Array& a, const <T&> b) {
  262.   <T>* u = a.x();
  263.   <T>* v = a.x() + a.n();
  264.   do {
  265.     while (u < v && *u $1 b)
  266.       u++;
  267.     if (u < v)
  268.       return 0;
  269.     u += a.l() - a.n();
  270.     } while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
  271.   return 1;
  272.   }
  273.  
  274. int operator $1 (const <T>Array& a, const <T>Array& b) {
  275.   if (a.n() == b.n()) {
  276.     if (a.m() == b.m()) {        // a.n() == b.n()
  277.       <T>* t = a.x();
  278.       <T>* u = b.x();
  279.       <T>* v = b.x() + b.n();
  280.       do {
  281.     while (u < v && *t $1 *u) {
  282.       t++;
  283.       u++;
  284.       };
  285.     if (u < v)
  286.       return 0;
  287.     t += a.l() - a.n();
  288.     u += b.l() - b.n();
  289.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  290.       return 1;
  291.       };
  292.     if (a.m() == 1) {            // a.n() == b.n() && a.m() != b.m()
  293.       <T>* u = b.x();
  294.       <T>* v = b.x() + b.n();
  295.       do {
  296.     <T>* t = a.x();
  297.     while (u < v && *t $1 *u) {
  298.       t++;
  299.       u++;
  300.       };
  301.     if (u < v)
  302.       return 0;
  303.     u += b.l() - b.n();
  304.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  305.       return 1;
  306.       };
  307.     if (b.m() == 1) {            // a.n() == b.n() && a.m() != b.m() && a.m() != 1
  308.       <T>* u = a.x();
  309.       <T>* v = a.x() + a.n();
  310.       do {
  311.     <T>* t = b.x();
  312.     while (u < v && *u $1 *t) {
  313.       t++;
  314.       u++;
  315.       };
  316.     if (u < v)
  317.       return 0;
  318.     t += b.l() - b.n();
  319.     u += a.l() - a.n();
  320.     } while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
  321.       return 1;
  322.       };
  323.     };
  324.   if (a.n() == 1) {            // b.n() != a.n()
  325.     if (a.m() == b.m()) {        // b.n() != a.n() == 1
  326.       <T>* t = a.x();
  327.       <T>* u = b.x();
  328.       <T>* v = b.x() + b.n();                        
  329.       do {
  330.     while (u < v && *t $1 *u)
  331.       u++;
  332.     if (u < v)
  333.       return 0;
  334.     t++;
  335.     u += b.l() - b.n();
  336.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  337.       return 1;
  338.       };
  339.     if (a.m() == 1) {            // b.n() != a.n() == 1 && a.m() != b.m()
  340.       <T>* u = b.x();
  341.       <T>* v = b.x() + b.n();
  342.       do {
  343.     while (u < v && *a.x() $1 *u)
  344.       u++;
  345.     if (u < v)
  346.       return 0;
  347.     u += b.l() - b.n();
  348.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  349.       return 1;
  350.       };
  351.     };
  352.   if (b.n() == 1) {            // a.n() != b.n() && a.n() != 1
  353.     if (a.m() == b.m()) {        // a.n() != b.n() == 1
  354.       <T>* t = b.x();
  355.       <T>* u = a.x();
  356.       <T>* v = a.x() + a.n();
  357.       do {
  358.     while (u < v && *u $1 *t)
  359.       u++;
  360.     if (u < v)
  361.       return 0;
  362.     t++;
  363.     u += a.l() - a.n();
  364.     } while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
  365.       return 1;
  366.       };
  367.     if (b.m() == 1) {            // a.n() != b.n() == 1 && a.m() != b.m()
  368.       <T>* u = a.x();
  369.       <T>* v = a.x() + a.n();
  370.       do {
  371.     while (u < v && *u $1 *b.x())
  372.       u++;
  373.     if (u < v)
  374.       return 0;
  375.     u += a.l() - a.n();
  376.     } while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
  377.       return 1;
  378.       };
  379.     };
  380.   a.error("nonconformant <T>Array $1 operands.");
  381.   return 0;
  382.   }
  383. ')dnl
  384. ifelse(basetype, complex, `',
  385. `int operator <  (const <T&> a, const <T>Array& b) {
  386.   <T>* u = b.x();
  387.   <T>* v = b.x() + b.n();
  388.   do {
  389.     while (u < v && a <  *u)
  390.       u++;
  391.     if (u < v)
  392.       return 0;
  393.     u += b.l() - b.n();
  394.     } while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
  395.   return 1;
  396.   }
  397.  
  398. comparison(`< ')
  399. ')dnl
  400. comparison(==)
  401. <T>Matrix operator & (const <T>Array& a, const <T>Array& b) {    // outer product
  402.   if (a.m() != b.m())
  403.     a.error("nonconformant <T>Array & operands.");
  404.   <T>* newx = new <T> [a.n()*b.n()];
  405.   <T>* p = newx;
  406.   <T>* s =  a.x();
  407.   <T>* t;
  408.   <T>* u;
  409.   <T>* v;
  410.   do {
  411.     t = b.x();
  412.     do {
  413.       u = s; v = t; *p = *u * *v;
  414.       while ((u += a.l()), (v += b.l()) < &b.x()[b.m()*b.l()])
  415.     *p += *u * *v;
  416.       p++;
  417.       } while (++t < &b.x()[b.n()]);
  418.     } while (++s < &a.x()[a.n()]);
  419.   return <T>Matrix(a.n(), b.n(), newx);
  420.   }
  421.  
  422. <T>Matrix operator ^ (const <T>Array& a, const <T>Array& b) {    // stack
  423.   if (a.n() != b.n())
  424.     a.error("nonconformant <T>Array ^ operands.");
  425.   <T>* newx = new <T> [(a.m()+b.m())*a.n()];
  426.   <T>* t = newx;
  427.   <T>* u =  a.x();
  428.   <T>* v =  a.x() + a.n();
  429.   do {
  430.     while (u < v)
  431.       *t++ = *u++;
  432.     u += a.l() - a.n();
  433.     } while((v += a.l()) <= &a.x()[a.m()*a.l()]);
  434.   u =  b.x();
  435.   v =  b.x() + b.n();
  436.   do {
  437.     while (u < v)
  438.       *t++ = *u++;
  439.     u += b.l() - b.n();
  440.     } while((v += b.l()) <= &b.x()[b.m()*b.l()]);
  441.   return <T>Matrix(a.m()+b.m(), a.n(), newx);
  442.   }
  443.  
  444. <T>Matrix operator | (const <T>Array& a, const <T>Array& b) {    // adjoin
  445.   if (a.m() != b.m())
  446.     a.error("nonconformant <T>Array | operands.");
  447.   <T>* newx = new <T> [a.m()*(a.n()+b.n())];
  448.   <T>* t = newx;
  449.   <T>* u =  a.x();
  450.   <T>* v =  a.x() + a.n();
  451.   do {
  452.     while (u < v)
  453.       *t++ = *u++;
  454.     t += b.n();
  455.     u += a.l() - a.n();
  456.     } while((v += a.l()) <= &a.x()[a.m()*a.l()]);
  457.   t = newx + a.n();
  458.   u =  b.x();
  459.   v =  b.x() + b.n();
  460.   do {
  461.     while (u < v)
  462.       *t++ = *u++;
  463.     u += b.l() - b.n();
  464.     t += a.n();
  465.     } while((v += b.l()) <= &b.x()[b.m()*b.l()]);
  466.   return <T>Matrix(a.m(), a.n()+b.n(), newx);
  467.   }
  468.  
  469. dnl 1$ = { ,*,/,+,-}
  470. define(assignment,
  471. `<T>Array& <T>Array::operator $1= (const <T&> b) {
  472.   <T>* s = X;
  473.   <T>* t = X + N;
  474.   do {
  475.     while (s < t)
  476.       *s++ $1= b;
  477.     s += L - N;
  478.     } while ((t += L) <= &X[M*L]);
  479.   return *this;
  480.   }
  481.  
  482. <T>Array& <T>Array::operator $1= (const <T>Array& a) {
  483.   if (N == a.N) {
  484.     if (M == a.M) {            // N == a.N
  485.       <T>* t = a.X;
  486.       <T>* u =   X;
  487.       <T>* v =   X + N;
  488.       do {
  489.     while (u < v)
  490.       *u++ $1= *t++;
  491.     t += a.L - a.N;
  492.     u +=   L -   N;
  493.     } while ((v += L) <= &X[M*L]);
  494.       return *this;
  495.       };
  496.     if (a.M == 1) {            // N == a.N && M != a.M
  497.       <T>* t;
  498.       <T>* u = X;
  499.       <T>* v = X + N;
  500.       do {
  501.     t = a.X;
  502.     while (u < v)
  503.       *u++ $1= *t++;
  504.     u += L - N;
  505.     } while ((v += L) <= &X[M*L]);
  506.       return *this;
  507.       };
  508.     };
  509.   if (a.N == 1) {            // N != a.N
  510.     if (M == a.M) {            // N != a.N == 1
  511.       <T>* t = a.X;
  512.       <T>* u =   X;
  513.       <T>* v =   X + N;
  514.       do {
  515.     while (u < v)
  516.       *u++ $1= *t;
  517.     t++;
  518.     u += L - N;
  519.     } while ((v += L) <= &X[M*L]);
  520.       return *this;
  521.       };
  522.     if (a.M == 1) {            // N != a.N == 1 && M != a.M
  523.       <T>* t = a.X;
  524.       <T>* u =   X;
  525.       <T>* v =   X + N;
  526.       do {
  527.     while (u < v)
  528.       *u++ $1= *t;
  529.     u += L - N;
  530.     } while ((v += L) <= &X[M*L]);
  531.       return *this;
  532.       };
  533.     };
  534.   error("nonconformant <T>Array $1= operands.");
  535.   return *this;
  536.   }
  537. ')dnl
  538. assignment(` ')
  539. assignment(*)
  540. assignment(/)
  541. assignment(+)
  542. assignment(-)
  543. <T>Array& <T>Array::operator >>= (int n) {
  544.   <T>* t = &X[M*L];
  545.   <T>* u = &X[M*L];
  546.   <T>* v = &X[M*L];
  547.   while (X <= (u -= L)) {
  548.     t -= L - N + n;
  549.     v -= L - N;
  550.     while (t > u)
  551.       *--v = *--t;
  552.     while (v > u)
  553.       *--v = *t;
  554.     };
  555.   return *this;
  556.   }
  557.  
  558. <T>Array& <T>Array::operator <<= (int n) {
  559.   <T>* t = X + n;
  560.   <T>* u = X;
  561.   <T>* v = X + N;
  562.   do {
  563.     while (t < v)
  564.       *u++ = *t++;
  565.     while (u < v)
  566.       *u++ = 0.0;
  567.     t += L - N + n;
  568.     u += L - N;
  569.     } while ((v += L) <= &X[M*L]);
  570.   return *this;
  571.   }
  572.  
  573. // Error Handling
  574. void default_<T>Array_error_handler(const char* msg) {
  575.   cerr << "Fatal <T>Array error. " << msg << "\n";
  576.   exit(1);
  577.   return;
  578.   }
  579.  
  580. one_arg_error_handler_t
  581.     <T>Array_error_handler = default_<T>Array_error_handler;
  582.  
  583. one_arg_error_handler_t set_<T>Array_error_handler(one_arg_error_handler_t f) {
  584.   one_arg_error_handler_t old = <T>Array_error_handler;
  585.   <T>Array_error_handler = f;
  586.   return old;
  587.   }
  588.  
  589. void <T>Array::error(const char* msg) const {
  590.   (*<T>Array_error_handler)(msg);
  591.   }
  592.  
  593. ifelse(basetype, complex, `', `
  594. <T>Matrix <T>Array::i(double epsilon) const {    // inverse
  595.   <T>*  a = new <T> [N*M];
  596.   <T>* uu = new <T> [M*N];
  597.   <T>* vv = new <T> [N*N];
  598.   <T>*  w = new <T> [N];
  599.   <T>** u = new <T>* [M];
  600.   <T>** v = new <T>* [N];
  601.   void svdcmp(<T>**, int, int, <T>*, <T>**);
  602.   for (int i = 0; i < M; i++)
  603.     u[i] = &(uu[N*i]);
  604.   for (int j = 0; j < N; j++)
  605.     v[j] = &(vv[N*j]);
  606.   for (i = 0; i < M; i++)
  607.     for (j = 0; j < N; j++)
  608.       u[i][j] = X[i*L+j];
  609.   svdcmp(u, M, N, w, v);    // Singular value decomposition.
  610.   <T> wmax = 0.0;        // Maximum singular value.
  611.   for (j = 0; j < N; j++)
  612.     if (w[j] > wmax)
  613.       wmax = w[j];
  614.   <T> wmin = wmax*epsilon;
  615.   for (int k = 0; k < N; k++)
  616.     if (w[k] < wmin)
  617.       w[k] = 0.0;
  618.     else
  619.       w[k] = 1.0/w[k];
  620.   for (i = 0; i < N; i++)
  621.     for (j = 0; j < M; j++) {
  622.       a[M*i+j] = 0.0;
  623.       for (k = 0; k < N; k++)
  624.     a[M*i+j] += v[i][k]*w[k]*u[j][k];
  625.       };
  626.  
  627.   delete [] w;
  628.   delete [] u;
  629.   delete [] v;
  630.   delete [] uu;
  631.   delete [] vv;
  632.  
  633.   return <T>Matrix(N, M, a);
  634.   }
  635. ')
  636. <T>Matrix <T>Array::t() const {        // transpose
  637.   <T>* newX = new <T> [N*M];
  638.   <T>* t = newX;
  639.   <T>* u;
  640.   <T>* v =    X;
  641.   do {
  642.     u = v;
  643.     do
  644.       *t++ = *u;
  645.     while ((u += L) < &X[M*L]);
  646.     } while (++v < &X[N]);
  647.   return <T>Matrix(N, M, newX);
  648.   }
  649.  
  650. <T>Matrix <T>Array::sum() const {
  651.   <T>* newX = new <T> [M];
  652.   <T>* t = newX;
  653.   <T>* u =    X;
  654.   <T>* v =    X + N;
  655.   do {
  656.     *t = *u;
  657.     while (++u < v)
  658.       *t += *u;
  659.     t++;
  660.     u += L - N;
  661.     } while ((v += L) <= &X[M*L]);
  662.   return <T>Matrix(1, M, newX);
  663.   }
  664.  
  665. <T>Matrix <T>Array::sumsq() const {
  666.   <T>* newX = new <T> [M];
  667.   <T>* t = newX;
  668.   <T>* u =    X;
  669.   <T>* v =    X + N;
  670.   do {
  671.     *t  = *u * *u;
  672.     while (++u < v)
  673.       *t += *u * *u;
  674.     t++;
  675.     u += L - N;
  676.     } while ((v += L) <= &X[M*L]);
  677.   return <T>Matrix(1, M, newX);
  678.   }
  679.  
  680. ifelse(basetype, complex, `
  681. #ifdef __ATT_complex__
  682. <T>Matrix <T>Array::map(const <T> (*f)(const <T>&)) const {
  683.   <T>* newX = new <T> [M*N];
  684.   <T>* t = newX;
  685.   <T>* u =  X;
  686.   <T>* v =  X + N;
  687.   do {
  688.     while (u < v)
  689.       *t++ = f(*u++);
  690.     u += L - N;
  691.     } while ((v += L) <= &X[M*L]);
  692.   return <T>Matrix(M, N, newX);
  693.   }
  694. #else
  695. <T>Matrix <T>Array::map(const <T> (*f)(      <T> )) const {
  696.   <T>* newX = new <T> [M*N];
  697.   <T>* t = newX;
  698.   <T>* u =  X;
  699.   <T>* v =  X + N;
  700.   do {
  701.     while (u < v)
  702.       *t++ = f(*u++);
  703.     u += L - N;
  704.     } while ((v += L) <= &X[M*L]);
  705.   return <T>Matrix(M, N, newX);
  706.   }
  707. #endif
  708. doubleMatrix <T>Array::map(const double (*f)(const <T>&)) const {
  709.   double* newX = new double [M*N];
  710.   double* t = newX;
  711.   <T>* u =  X;
  712.   <T>* v =  X + N;
  713.   do {
  714.     while (u < v)
  715.       *t++ = f(*u++);
  716.     u += L - N;
  717.     } while ((v += L) <= &X[M*L]);
  718.   return doubleMatrix(M, N, newX);
  719.   }
  720. #ifndef __ATT_complex__
  721. doubleMatrix <T>Array::map(const double (*f)(      <T> )) const {
  722.   double* newX = new double [M*N];
  723.   double* t = newX;
  724.   <T>* u =  X;
  725.   <T>* v =  X + N;
  726.   do {
  727.     while (u < v)
  728.       *t++ = f(*u++);
  729.     u += L - N;
  730.     } while ((v += L) <= &X[M*L]);
  731.   return doubleMatrix(M, N, newX);
  732.   }
  733. #endif
  734. ',
  735. `<T>Matrix <T>Array::map(const <T> (*f)(      <T> )) const {
  736.   <T>* newX = new <T> [M*N];
  737.   <T>* t = newX;
  738.   <T>* u =  X;
  739.   <T>* v =  X + N;
  740.   do {
  741.     while (u < v)
  742.       *t++ = f(*u++);
  743.     u += L - N;
  744.     } while ((v += L) <= &X[M*L]);
  745.   return <T>Matrix(M, N, newX);
  746.   }
  747.  
  748. <T>Matrix <T>Array::min() const {
  749.   <T>* newX = new <T> [M];
  750.   <T>* t = newX;
  751.   <T>* u =    X;
  752.   <T>* v =    X + N;
  753.   do {
  754.     *t = *u;
  755.     while (++u < v)
  756.       if (*t > *u)
  757.     *t = *u;
  758.       t++;
  759.       u += L - N;
  760.     } while ((v += L) <= &X[M*L]);
  761.   return <T>Matrix(1, M, newX);
  762.   }
  763.  
  764. <T>Matrix <T>Array::max() const {
  765.   <T>* newX = new <T> [M];
  766.   <T>* t = newX;
  767.   <T>* u =    X;
  768.   <T>* v =    X + N;
  769.   do {
  770.     *t = *u;
  771.     while (++u < v)
  772.       if (*t < *u)
  773.     *t = *u;
  774.     t++;
  775.     u += L - N;
  776.     } while ((v += L) <= &X[M*L]);
  777.   return <T>Matrix(1, M, newX);
  778.   }
  779.  
  780. int <T>Array::min_index() const {
  781.   <T>* t = X;
  782.   <T>* u = X;
  783.   <T>* v = X + N;
  784.   do {
  785.     do
  786.       if (*t > *u)
  787.     t = u;
  788.       while (++u < v);
  789.     u += L - N;
  790.     } while ((v += L) <= &X[M*L]);
  791.   return (t - X);
  792.   }
  793.  
  794. int <T>Array::max_index() const {
  795.   <T>* t = X;
  796.   <T>* u = X;
  797.   <T>* v = X + N;
  798.   do {
  799.     do
  800.       if (*t < *u)
  801.     t = u;
  802.       while (++u < v);
  803.     u += L - N;
  804.     } while ((v += L) <= &X[M*L]);
  805.   return (t - X);
  806.   }
  807. ')dnl
  808. ifelse(basetype, complex, `
  809. <T>Matrix fft(const <T>Array& a) {
  810.   extern void four1(double*, int, int);
  811.   <T>Matrix    result(a);
  812.   for (int i = 0; i < result.m(); i++)
  813.     four1((double*) result[i] - 1, result.n(), 1);
  814.   return result;
  815.   }
  816.  
  817. <T>Matrix fft(const doubleArray& a) {
  818.   extern void realft(double*, int, int);
  819.   <T>Matrix    data(a.m(), a.n());
  820.   doubleArray(2*a.n(), a.m(), a.n(), (double*) data.x()) = a;
  821.   for (int i = 0; i < a.m(); i++) {
  822.     realft((double*) data[i] - 1, a.n()/2, 1);
  823.     for (int j = 1; j < a.n()/2; j++)
  824.       data[i][a.n()-j]    =  conj(data[i][j]);
  825.     data[i][a.n()/2]    =  imag(data[i][0]);
  826.     data[i][0]        =  real(data[i][0]);
  827.     };
  828.   return data;
  829.   }
  830. ', `
  831. <T>Matrix ffct(const <T>Array& a) {
  832.   extern void cosft(double*, int, int);
  833.   <T>Matrix    y(a);
  834.   for (int i = 0; i < a.m(); i++)
  835.     cosft(y[i] - 1, a.n(), 1);
  836.   return y;
  837.   }
  838. ')
  839. // Constructors
  840. <T>Matrix::<T>Matrix(const <T>Matrix& a) :
  841. <T>Array(a.M, a.N, new <T> [a.M*a.N]) {
  842.   <T>* t = a.X;
  843.   <T>* u =   X;
  844.   <T>* v =   X;
  845.   while ((v += N) <= &X[M*N]) {
  846.     while (u < v)
  847.       *u++ = *t++;
  848.     t += a.L - a.N;
  849.     };
  850.   }
  851.  
  852. <T>Matrix::<T>Matrix(const <T>Array& a) :
  853. <T>Array(a.M, a.N, new <T> [a.M*a.N]) {
  854.   <T>* t = a.X;
  855.   <T>* u =   X;
  856.   <T>* v =   X;
  857.   while ((v += N) <= &X[M*N]) {
  858.     while (u < v)
  859.       *u++ = *t++;
  860.     t += a.L - a.N;
  861.     };
  862.   }
  863.  
  864. ifelse(basetype, complex, `
  865. <T>Matrix::<T>Matrix(const doubleArray& r, const double i) :
  866. <T>Array(r.M, r.N, new <T> [r.M*r.N]) {
  867.   <T>* u =   X;
  868.   <T>* v =   X;
  869.   double* t = r.X;
  870.   while ((v += N) <= &X[M*N]) {
  871.     while (u < v)
  872.       *u++ = complex(*t++, i);
  873.     t += r.L - r.N;
  874.     };
  875.   }
  876.  
  877. <T>Matrix::<T>Matrix(const doubleArray& r, const doubleArray& i) :
  878. <T>Array(r.M, r.N, new <T> [r.M*r.N]) {
  879.   if(i.M == M && i.N == N) {
  880.     <T>* u = X;
  881.     <T>* v = X;
  882.     double* t =  r.X;
  883.     double* k =  i.X;
  884.     while ((v += N) <= &X[M*N]) {
  885.       while (u < v)
  886.     *u++ = complex(*t++, *k++);
  887.       t += r.L - r.N;
  888.       k += i.L - i.N;
  889.       };
  890.     }
  891.   else
  892.     error("nonconformant <T>Array operands.");
  893.   return;
  894.   }
  895. ')dnl
  896. <T>Matrix::~<T>Matrix()
  897. { delete [] X; }
  898. popdef(`index')dnl
  899. popdef(`shift')dnl
  900. popdef(`format')dnl
  901.  
  902.